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ABSTRACT 


The transient surface temperature distribution is 
determined for the flat plate and sphere subjected to cooling 
by combined convection and radiation. In the study, the 
initial boundary value problem is reduced to a singular 
nonlinear Volterra integral equation of the second kind using 
the integral transform method. Several numerical techniques 
are introduced in an attempt to find an approximate solution 
of the problem: The method of successive approximations, the 
Runge-Kutta method, and the finite difference method. The 
integral equation is solved numerically by the Runge-Kutta 
method of orders 1, 3, and 5. In addition, the finite 
difference method is implemented to solve the initial boundary 
value problem, and the solutions are seiewied with those 
generated by the Runge-Kutta method. All the numerical 
results are presented graphically. Limitations and 
difficulties involved in these schemes are discussed. At the 
end, a numerical algorithm for solving the problem is 


proposed. 
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I. ANALYTIC SOLUTIONS OF THE HEAT EQUATION SUBJECT TO 


CONVECTIVE AND RADIATIVE BOUNDARY CONDITIONS 


A. INTRODUCTION 

During the 60's, space technology advanced so much that 
the research of the temperature behaviour of bodies exposed to 
a deep space environment became crucial. In particular, 
transient heat or cooling of solids of different shapes by 
convection and thermal radiation was becoming highly important 
in many engineering applications. An example of these 
applications is the temperature distributions of rocket 
motors. An extensive investigation of the problem has been 
conducted and a lot of literature on the subject was published 
during the 60's and 70's. A detailed review of most of these 
papers is not intended here; instead a brief summary of the 
major ones will be given. 

As early as 1962, Fairall, et.al.({6] generated a numerical 
solution for the problem using an explicit finite difference 
scheme; this paper served as pioneer work in the area of the 
research. Later, various finite difference schemes were 
devised to deal with the nonlinear boundary condition. The 
main difficulty in these schemes is the appearance of severe 
oscillations in the determined temperature values for high 


heat flux situations. Von Rosenberg [10] proposed a hybrid of 


an iterative technique and implicit finite difference schemes 
to deal with the nonlinear boundary condition. On the other 
hand, Crosbie and Viskanta [3,4] transformed the governing 
equations into a nonlinear Volterra integral equation of the 
second kind and applied the method of § successive 
approximations to solve the integral equation. Milton and 
Goss [8,9] developed some heuristic stability criteria for 
explicit finite difference schemes with nonlinear boundary 
conditions. It turns out that a very restrictive time step is 
required for numerical stability which may result in requiring 
a prohibitive amount of computer time to calculate the long 
time evolution of the solutions. Williams and Curry [12] 
surveyed several methods for treating the nonlinear boundary 
condition in implicit schemes and compared their accuracy and 
efficiency. 

Nonlinearity is commonplace in natural phenomena. 
Unfortunately, a nonlinear problem often doesn't lend itself 
to a closed form solution. The problem of transient heat- 
conduction in a solid becomes nonlinear when the surface of 
the body is subjected to thermal radiation. When energy 
transfers through the wall of a body, two cases arise: 
convection and thermal radiation. The convective heat 
transfer describes the situation where heat is dissipated 
according to Newton's Law of Cooling, which states that the 
rate at which heat is transferred from the body to a 


surrounding is proportional to the difference in temperature 


between the body and the environment. The boundary condition 
that describes convection is nonlinear except for the case 
where the heat-transfer coefficient is independent of surface 
temperature, which is technically called forced convection. 
The radiative heat transfer is based on the Stefan-Boltzmann 
Law, which states that the heat flux is proportional to the 
difference between the surface temperature to the fourth 
power and the source temperature. Pure radiation or pure 
convection occur whenever one mode of energy transfer 
predominates over the other. 

It is the purpose of this thesis to consider the one- 
dimensional transient heat conduction problem resulting from 
a combined convective and radiative heat flux with the 
objective of determining the surface temperature fields using 
the numerical methods which are discussed in this study. 
Another purpose of this thesis is to explore the limitations 
and difficulties involved in these schemes. References to the 
work done in similar areas are presented to allow the reader 
further investigation. 

Analytic solutions are derived in one dimension. However, 
the resulting solutions are not in closed form, and thus 
impractical to use. Hence, numerical techniques will be 
studied and employed in the computer in an attempt to find an 
approximate solution. Numerical results, found by 
implementing some of the numerical methods discussed below, 


will be presented and compared. In the conclusion, a 


numerical scheme is proposed as an alternative to the existing 
methods. It is open to the readers for justification. 

Sections 1(C) and 1(D) describe the derivation of the 
integral representations of the one dimensional transient heat 
conduction problem subjected to a combined convective and 
radiative boundary condition in a rectangular coordinate 
system. Two integral transform methods, namely the Laplace 
transform and the eigenvalue expansion, are presented. 
Observation and comparison are made for the integral equations 
to yield some useful information about the solutions. 

In Chapters II and III, numerical methods for the 
solutions of the nonlinear Volterra integral equations of the 
second kind are described. In particular, the method of 
successive approximations and the Runge-Kutta method are 
outlined in detail. A brief remark is given for their 
advantages and limitations in finding solutions to the 
integral equation. 

Chapter IV describes a numerical method which is directly 
applied to the governing partial differential equation. The 
technique is called the finite difference method. It is 
basically a hybrid of finite difference techniques and an 
iterative scheme proposed. A suggestion is made for the 
improvement of the algorithn. 

In Chapter V, numerical results produced by some of the 
discussed numerical schemes are presented. The implementation 


of various methods gave a practical sense of their advantages 


and limitations. Graphs and tables are set up in such a way 
that a comparison can be made. 

In the next section, a statement of the problem is given. 
In the statement, the basic assumptions, the governing 


equation and the boundary-initial conditions are included. 


B. STATEMENT OF THE PROBLEM FOR OBTAINING THE SURFACE 
TEMPERATURE 
Considering the one-dimensional, transient, conduction 
heat transfer problem with combined convection and radiation 
at its surface, the following assumptions have been made: 


1. One-dimensional heat transfer to a solid of a finite 
length. 


2. The solid medium is pure, isotropic, homogeneous, and 
Opaque to thermal radiation. 


3. All thermodynamic and transport properties are 
independent of temperature. 


4. The solid does not contain any heat sources or sinks. 
5. The fluid is transparent to thermal radiation. 


6. The fluid temperature and the ambient temperature are 
constant. 


The non-dimensional form of the governing partial 


differential equation for the temperature U(x,t) and the 


appropriate initial boundary conditions are 


ORE = O<x<1, t>0; (1.1) 


with initial condition 


U(x,0) = g(x) (1.24) 

and boundary conditions 
(i) a iiall(® -f - @, U(0,t) = 0 (1.2b) 
(ii) ka ras - @,U(1,t) = -hU‘(1,t). (1.2c) 


Note: @, and @, can be any real number, except both cannot be 
zero at the same time. «a; 1S a non-zero real number, and h is 
a positive real number. 

The next section will deal with solving the partial 
differential equation (1.1) with initial and boundary 
conditions (1.2a-c) by the Laplace transform method and the 
eigenvalue expansion method. As an illustration, two special 
cases with specific values of @,, @, @3, and h will be 
considered, and the analytic solutions of these cases at the 
surface will be derived. It will be shown that the surface 
temperature satisfies a singular Volterra integral equation of 
the second kind. At the end of the chapter, we will present 
the solutions and indicate some useful information about the 


integral equations. 


C. THE LAPLACE TRANSFORM METHOD 

In this section, the Laplace transform of equation (1.1) 
with associated boundary conditions (1.2b,c) is first obtained 
with respect to time. The resulting boundary value problem is 
in terms of the Laplace transform of the required solution. 
Next, the equations are solved for the transformed 
temperature, and the solution of the stated problem can be 
found by taking the inverse Laplace transform of the 
transformed solution. From experience, it can be expected, 
the Laplace inversion is of some difficulty. To simplify the 
situation, specific values of «@,, &@,, &@,, and h are considered 
so that the inverse process is practical without loss of 
generality. It should be noted that there does exist an 
inverse Laplace transform for other cases of a more general 
nature. 

Now, define the transform of the temperature function, 


U(x,t), with respect to time as follows 
or u(x, t))(s) = [-a(x, ty estat = U(x,S). (1.3) 


After the transformation, the temperature function becomes a 
function not only of x but also of the parameter s. Assuming 
that the derivatives with respect to x pass through the 
transform (differentiation can be accomplished before 


integration), we have 


st Slee = a = ae (1.4) 


07u(x,t) MR OU CO pey. -st _ Oxy 
SE mS sees dt a (1.5) 


The rule for transforming a derivative with respect to time 
can be found using integration by parts. Thus, the Laplace 
transform of the derivatives of U(x,t) with respect to the 


transformed variable t is given by 


a eee ie = jg a: = sSU(x,S) - U(x,0). (1.6) 


Now, applying the Laplace transform to the initial-boundary 
value problem (1.1), (1.2a-c) we remove all time derivatives. 
Holding s fixed, we have the following ordinary differential 


equation in x 


d*U(x,S) 


a2 - SU(x,S) = -g(X), 0<x<1 (1.7) 


with boundary conditions 


- @,U(0,S) =0, for x =0 (1.84) 


. du(o, s) 


= ~ @,U(1,S) = -h&[U‘(1,s)], for x=1 (1.8b) 


Notice that the initial condition, g(x), is incorporated in 
the ordinary differential equation. In order to solve (1.7) 
and (1.8a,b), we must first solve for the general solution of 
the corresponding homogeneous differential equation and a 
particular solution of (1.7) satisfying (1.8a,b). Now, 
consider the general solution of the homogeneous equation for 


ae7), 
U,o.n(X,S) = Ae** + Be**, (1.9) 


where 


4,, =+y7s (1.10) 


W-s=0. (1.11) 


In the following paragraph we employ the method of variation 


of parameters to solve for a particular solution of (1.7). 


Let 
U,(X,S) = ULV, Case v(x, Ss); (1.12) 


be a particular solution where U,(x,s) and U,(x,s) are any two 
linearly independent solutions oof the corresponding 
homogeneous equation. In this case, choose U,(x,s) = e’* and 
U,(x,S) = e’*, The object here is to find v,(x,s) and v,(x,s) 


such that the following equations are satisfied 
evaxy’! (x,s) + e7 Bx! (x,s) = 0, (1.13) 


vse’? *vi(x,s) - ¥seV? *vi(x,s) =- g(x) (1.14) 


By Cramer's rule, 


vi (x, S) = Sieiein= (1.15) 
“28 
and 
Vi 
v(x, s) = gixlev™ | (1.16) 
2/s 
By integrating (1.15) and (1.16), we obtain 
-V5 
v, (X, S) - -[*Gizle wa + v, (0,8) (ttl: 7 >) 
0 2/s 
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_ f*g(z)evdz _ 
v5 (x, Ss) —as v, (0,3) (1.18) 


Thus, the general solution to (1.7) and (1.8a,b) is 


U(X, S) = U,(X,S) + U,(X,S) (1.19) 


that is, 


U(x,s) = Aev® + Be-V8 + eve, (x,s) + e-V8*v, (x, 5), (1210) 


where A and B are arbitrary constants and vU,(X,S), UVU,(x,sS) are 
given by (1.17) and (1.18), respectively. To determine A and 
B, boundary conditions (1.8a,b) are used along with ene 
following procedure. The derivative of U(x,s) from equation 


(1.20) is found to be 


Stes) = A/sev™ - B/se-v=* + Vsev®yv, (x, s) + 
every) (x,s) + e7 Bxy! (x,s) - Jsevv,(x,S). (1.21) 


Let x = 0. Then (1.15) and (1.16) give 


vi(0,s) +v5(0,s) =0. 


(1.18a), (1.20), and (1.21) then imply 


11 


a, [AVS - BYS + Y5v,(0,s) - /Sv,(0,5)] 


- @,[A + B +v,(0,s) + v,(0,sS)] = 0 (1 Ze 


By rearranging the terms, (1.22) becomes 
A(a@.jS - @,) - Bla./s + a) = 


(a,J/S + @)v,(0,S) - (a,/S - a@)v,(0,5S). (1223) 


Similarly let x = 1. Then (1.15) and (1.16) give 


evéy'(1,s) + evVavi(1,s) =0. (1.24) 


Therefore (1.8a), (1.20), and (1.21) then imply 
A/sev® - B/se-v5 + /sevV8v,(1,s) - /se-v5v,(1,8s) - 
a,[Aev® + Be-v8 + evV8v, (1,5) + e-V8v,(1,s)] = 


—Weae (1, Coe (1.239 


By a similar manipulation of the terms, (1.25) becomes 
A(/sev® - a,ev*) - B(/se-v* + a,e-V%) = | 
(/se-v* + a,e-¥%)v, (1,5) - (¥Sev® - a,e¥*)v,(1,s) - 


Hee (dat)de (1.26) 


12 





Equations (1.23) and (1.26) form a system of two equations in 
the two unknowns A and B. By Cramer's rule, A and B are as 


follows 


num1 
aa ee) 
den ( ) 





where 
num1 = 
{(a@,/S-a,) v, (0,5) +(a,/S+a,) [v, (1, 5) -v, (0,8) ]}(Vse-vFt+a,e-V*) 


a [(a, St+a,)v,(1,5s)] (/sev*-a,ev*) = hL[ut(1, b) ] (a, Sta.) ’ 





and 
den = 
(a, - @,0,)/s(ev® + eV?) + (a,s - «@,0a,) (eV8-e-v*) 
=p Tn (Gos) 
den 
where 


num2 = 
{(a,/S-a,) [v,(0,s) -v, (1,5) ]-(a,/st+a,) v, (0,5) }(¥sev5-a,ev8) 


- HL [U*(1, t)] (ayVS-a,) + [ (a,V5-a) v2 (1,8) ] (VSeVF+a,€-¥9) 


and den is the same. 
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Thus, the general solution of (1.7), (1.8a,b) is given by 
(1.20) where A, B, U,(xX,S), and vU,(x,S) are given by (1.27), 
(1.28), (1.17), and (1.18), respectively. Theoretically, the 
analytic solution of partial differential equation (1.1) with 
initial and boundary conditions (1.2a,b,c) can be obtained by 
taking the Laplace inversion of U(x,s), and thus, the surface 
solution can be found by putting x equal to one in U(x,s). In 
practice, however, the inverse Laplace transformation process 
is highly unstable in that singularities may exist. Also, the 
transforms are difficult to find. In the next paragraph below 
we consider two special cases where the inversion is feasible. 
In each case, values for the parameters correspond to a 


specific geometrical configuration of a body. 


Case 1: @,=1, @ = 0, @ =-1, h=1 

This- set of values corresponds to a "flat plate" with a 
given initial temperature and which is being heated or cooled 
by combined convection and radiation. The term "flat plate" 
is taken here to mean a solid slab of finite thickness which 
is bounded by a pair of vertical lines at + 4% thus having a 
width of 1. Substituting the given values for the parameters 


in (1.20) we obtain the transformed surface temperature 


U(1,s) = Aev® + Be-V8 + ev8y, (1,5) + e-¥8v,(1,8), (1.29) 


where, 


14 


/s(ev8 + e-v8) + s(ev4 J e-v3) 


ieemeenL ees = oF) neo" (1. ¢) lve 


(laad'O ) 
Vs(ev® + e-¥8) + s(ev® - e-v®) 


(4. - a,0,)¥s(ev* + ev") + (a,5 - a3a,) (ev* - ev") 


, 7 belt (2, t)1 V5 + [YSv2(1,s)] (¥Se%? - eV) 


(1.31) 
Kes m a,0,) /s(evs oa e-v8) + (a8 = 0505) (evs 5 @-v5) 


_.. f.gix'he-Ve 
v, (1,8) ti a dx’ + v, (0,8) (ane '2 ) 


and 


Par) = jf, cien ax! + v,(0,5) 


(1.33) 
2/s 


Now substituting (1.30)-(1.33) into (1.29) and simplifying the 


results gives: 


Rg (x!) eV®/dx! + fig (x!) e-v#'dx! 
(evs + e-vé) +} Vs (ev4 a e-v8) 


U(1,s) = 


pees Wet, EN (eve + eV) _ 


(1.34) 
(ev? + e-v4) + JS(eV8 - e-V8) 
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Suppose the initial temperature is 1, that is, 


g(x) = 1. (1.35) 


The boundary conditions associated with the given values of 
the parameters and initial condition (1.35) constitute a 
cooling process. With (1.35), the transformed surface 


temperature becomes 


1 (ev? - e-v4) - h f[u4(1, t)] (ev® + e-v4) 


os) ae . 4.38) 
(ev? + e-v®) + /s(e/s - e-V*) 


If (1.36) is multiplied through by 


_1 “) SY tage’ ot /SNERE NE) 


3 (ev - e-v8) 


and then simplified, 


S[h U4(1, t) + U(1, t)] (ev® + e-v®) 


(1.37) 
VS(ev® - e-v8) 


U(1l,s) = - 
Ss 


is obtained. Equation (1.37) is ready to be inverted. In 
order to perform the inversion of (1.37), the following two 
Laplace transforms have to be computed 


(eR + ems) 


g-2( 2) and | 
= VE Claas 
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In fact the transforms can be found from any standard Laplace 
transform table. By the convolution theorem, the surface 


temperature in time t is given by 


ck* 
igi + 2 is (t=t) 
u(1,t) =1 =f Rens OT ota) + U(1,t)]dt (1.38' 


0 vere — t) 


and by the Poisson summation formula [14], (1.38) can be 


written as 


U(1,t) =1 or GR MO Saldaaels [u*(1,t)+U(1,t)] dt. (1.39) 


Hence, the problem of transient cooling of a flat plate by 
combined convection and thermal radiation has been reduced to 
solving a nonlinear Volterra integral equation of the second 


kind. - 


Case 2: 4, = 0, @ = rl, @4;=1, h=1 

This set of values corresponds to the case where a 
spherical body of radius 1 with a given initial temperature 
is being heated or cooled by combined convection and 
radiation. Since the procedures used to solve the problem are 
basically those described in case 1, the mathematical details 
will be omitted and only the main steps will be presented. 
Consider equation (1.20), the general solution of the boundary 


value problem. The given values for the parameters are first 
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substituted into (1.17), (1.18), (1.27), and (1.28). Then, 
(1.20) is simplified as in the previous case. After a tedious 


calculation, the transformed surface temperature is given by 


footx e- Bx! xe! = foot) ev Ex! x! + 
ms 0 


U(1,8s) 
(ev? - e-v®) - /s(ev? + ev) 


h £[u4(1,t)] (ev? + e-v®) 


+ ) (1.40) 
(ev® - e-v8) - Ys(ev® + e-v®) 
Suppose the initial temperature is chosen to be 
g(x) = xX. (1.41) 


Boundary conditions associated with the given values of the 
parameters and initial condition (1.41) again constitute a 
cooling process. With (1.41), the transformed steeds 
temperature becomes 


h @[(u4(1, t)] (ev® - e-v8) (aap) 


1 
U Le = — = 
( s) S /s(ev8 ne e-v8) a (evs = e-vé) 


which is now ready to be inverted. In order to perform the 
Laplace inversion of (1.42) the following two inverse Laplace 
transforms need to be computed 


(ev® - ev) 


eoet fo] and ot [ —————___-$ = —__— 
dl /S(ev3 + e-v8) - (ev? - e-v4) 
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The first inverse Laplace transform is obvious. However, the 
second one is not so obvious. Details of the derivation of 
the second inverse Laplace transform are given in [1]. The 
surface temperature in time t, obtained by inverting (1.42), 


is 


D(1,£) = 1 - f°[3 + 20, eM] uta, 2) ae, (1.43) 


where f, is the k® positive root of the transcendental 
equation 

B, = tan f, . (1.44) 
Hence, the problem of transient cooling of a sphere by 
combined convection and thermal radiation has been again 
reduced to solving a nonlinear Volterra integral equation of 
the second kind. As we have mentioned above, one of the 
Grawbacks of the Laplace transform method is that there are 
only a few cases in which the transformed solution can be 
practically inverted into the required solution. In the next 
section, the eigenvalue expansion method is introduced as an 
alternative to the above method. One may find the eigenvalue 
method more practical for solving for the analytic solution of 


the heat equation with nonlinear boundary conditions. 


D. THE EIGENVALUE EXPANSION METHOD 
The fundamental idea of the eigenvalue expansion method is 


to transform the given boundary value problem by the 
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eigenfunctions obtained from the associated eigenfunction 
problem. By the completeness theorem (which states that any 
piecewise smooth function can be represented by a generalized 
series of eigenfunctions) we can show that separation of 
variables, i.e., u(x,t) = X(x)T(t), may lead to the solution 
of the problem expressed as an infinite sum of the 
eigenfunctions with appropriate coefficients determined by the 
orthogonality property of eigenfunctions. Applying these 
procedures to the partial differential equation (1.1) and 
initial boundary conditions (1.2a,b,c) yields the following 


main results 


+ B2X(x) = 0, O0<x<1 (1.45) 


with boundary conditions 


= 





a, SEO) - @,X(0) = 0, (1.46) 
and 
eee ~@,X(1) =0. (1.47) 


Parameters @, and @, can be any real number except they cannot 
be zero at the same time. @,; 1S a non-zero real number. 
According to the theory of ordinary differential equations, 


the general solution of (1.45) is 
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X(x) = c,cos(Bx) + c,sin(Bx) . (1.48) 


Applying boundary conditions (1.46) and (1.47) to equation 


(1.48) gives the following system of equations 


“Bc, = a,c, (1.49) 


(ep ea, cesp = (Cc/pe+ c.g,)sinB™ (2rso) 


Note that boundary value problem (1.45 - 1.47) is in the 
class of Sturm-Liouville problems for which all eigenvalues 
are real and the eigenfunctions corresponding to different 
eigenvalues are orthogonal. Thus, if the parameters in (1.49) 
and (1.50) are specified, there will exist eigenvalues, 6,, 
where n = 1,2,..., and the corresponding eigenfunctions, 
X,(x), such that the temperature function, U(x,t), can be 


expanded in a Fourier expansion of the form 


U(xXptho= Woo Uplt) X,{x), (dy. Side) 


where the Fourier coefficients, U,(t), are given by 


U,(t) = [roe t) X,(x) dx . (1.52) 
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Now, taking the finite Fourier integral transform of the heat 
equation (1.1) with respect to X,(x) gives 


di _ £1 &U 
“ef U(x t) X_(x) dx = ff Sa ka(*) ax . (1.53) 


Performing integration by parts of the right hand expression 
in equation (1.53) and substituting in (1.52) yields the 
following ordinary differential equation for U,(t) 


du, (t) a 0oU(1,t) 


_ OU{(0, t) - / 
= a Ky (1) - ASX, (0) - U(1, €) Xn (1) 


+ U(0, t) X4(0) + [-Ulx, t) x(x) dx . (1.54) 


With boundary conditions (1.46) and (1.47), the right hand 
side of (1.54) can be simplified. Then, by the integrating 
factor method, the solution of equation (1.54) can be 
obtained. Hence, the resulting integral equation for U(x,t) 
takes the form of (1.51) with U,(t) solved in (1.54). Lastly, 
by putting x = 1, a nonlinear Volterra integral equation of 
the second kind for the surface temperature U(1,t) is 
obtained. 

As in the previous section, the integral equation for the 
surface temperature will be explicitly determined for two 
special cases: the flat plate and the sphere. Details of the 


derivation of the solution will be produced in the case of the 


22 


flat plate, but only major results will be given in the case 
of the sphere. 
Case 1: @, =1, @ =0, @, = -1, h=1 

As mentioned in section 1(C), this set of parameters 
corresponds to the geometrical configuration of a flat plate. 
Substituting the values of a,, @,, and a, in (1.49) and (1.50), 


Cc, equals zero, and (1.50) leads to 
cosB, = B,sinB, => om = tanB, , (gae5 5) 
n 
where cos f, # 0. 
So, the family of orthogonal eigenfunctions are 


X,(x) = cos(B,x) , (1.56) 


where n = 1,2,3,..., and {f,},..° is the set of distinct 


eigenvalues which are the roots of (1.55) with the property 


0<B,<B,<B,<... 


Next, applying the finite Fourier integral transform of the 
heat equation (1.1) yields (1.54) in terms of X,(x). Using 


the boundary conditions 


Wnt aa 


ax ’ ( te, 57) 
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0U(1,t) 


o se Se a) ie I ae la a I el (1.58) 
X"(1) aaa = 07. (1.59) 
xa), = Ome (1.60) 
and the fact that 
Xi(x) = — Pex (=) (1.61) 


produces the following ordinary differential equation for 
U, (t) 


ae + B2U,(t) = - bX,(1) U4 (1p) 


(1.62) 


Note that (1.62) is a first order linear ordinary differential 


equation. We find the solution to be 


u(t) = U,(0)e Pat - b x, (1) feet 8 (2, 2) dt, (1.63) 


where 


y_(0) = fg x) X,(x) dx . (1.64) 
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Thus, with h = 1, the integral equation for U(x,t) takes the 


form 


[u,(0) e-P#® - (1) x, (1) fe PH" Ut (2,2) dt] x5 (2) 


[x (x) dx 


where U,(0) and X,(x) are defined by (1.64) and (1.56), 
respectively. Lastly, by putting x = 1, the integral solution 


for the surface temperature U(1,t) is determined to be 


e Paty (1) fog X,(x) dx 


Aa, 


ieee) = ae 


n=l 


[-% (x) dx 


=f weit) eo (2, 1) ae 
2 
[, x (x) dx 


where g(x) is the initial condition, and X,(x), and 6, are 


defined as above. 
Case 2: @, = 0, @=-r1, @&= 1, hH=l 
In this case, a spherical body is considered. In a 


Similar fashion, the family of orthogonal eigenfunctions can 


be found and are given by 


X,(x) = sin(B,x) , (1.66) 
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where n= 1,2,3,..., and f§, is the set of distinct eigenvalues 


that are the roots of 


B, = tan 6, (1.67) 


with the property 


O < 8. <eb) <5b2< 22. : 


After applying the finite Fourier integral transform of heat 
equation (1.1) with respect to X,(x), the following ordinary 


differential equation for U,(t) is obtained 


du,(t) 


ae B2u5(t) = -et, (a)o* Gee (1.68) 


Thus, the solution of equation (1.68) is 


u.(t) = u,(0)e Ft - b x, (1) [oe P#**) 4 (1,2) dt, (1.69) 


where 


U,(0) = [ g(x) X,(x) dx . (1.70) 


So therefore, with h = 1, the integral equation for U(x,t) 
takes the form 


[y,(0)eP* - (2) x, (2) fe PP ot (a,c) ae) xq (x) 
B(x, t) = 0 ($e 


n=1 


fuxe (x) dx 
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where U,(0) and X,(x) are defined by (1.70) and (1.66), 
respectively. Lastly, by putting x = 1, the integral equation 
for the surface temperature U(1,t) becomes 


_ 2 Fx, (1) f g(x) x, (x) ax 
u(1,t) =  {——___-* ________ 


n=1 


[ Xa (x) dx 


= [7a (2) oF 4 (1, 2) dt 
0 


[x (x) dx 


where g(x) is the initial condition, and X,(x) and §, are 


defined as above. 


E. REMARKS 

The solution presented above is not complete in the sense 
that the surface temperature is only determined for two cases. 
The solution for other geometrical configurations can be found 
in some of the literature listed in the references, 
specifically 3, 5, 6, and 11. 


The surface temperature solutions which have been derived 


above by both methods fall into the form 


U(1,t) = o(t) - Af la+ Dy de FM P[O(1, 2) de, (1.72) 


where F is a nonlinear function of U(1,t), and c,, b,, a, and 


h are some constants. Equation (1.72) is a nonlinear 
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Volterra integral equation of the second kind. 9$(t) is a 
function which is usually called the "lag" part of the 
integral equation. The integral in (1.72) is often referred 
to as the "Volterra" part of the integral equation. In 
addition, the piece within the braces of the Volterra part is 
called the "kernel" of the integral equation. As these 
integral equations are being examined, several facts about 
(1.72) are summarized as follows: 
1). All of these integral equations are singular because as 
tT approaches t, the kernel blows up to infinity. 
2%: All of the infinite series satisfy the following 
property: 

If f(t-t) is used to denote an infinite series, 
then lim... f(t) = constant, thus remaining finite. 
3) The lag part, o(t), and the kernel of the integral 
equations are determined by the geometry of the body 
considered. 

The above "facts" are concluded from the two special cases 
without loss of generality. In each of the next three 
chapters, a different numerical method for solving the problem 
stated in section A will be introduced. Both the method of 
successive approximations and the Runge-Kutta method are 
numerical techniques used to deal with the integral 
representation of the problem, whereas the finite difference 


method is applied directly to the governing equations. 
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II. THE METHOD OF SUCCESSIVE APPROXIMATIONS 


A. INTRODUCTION 

The surface temperature of a body subject to a combined 
convective and radiative boundary condition, as seen in 
Chapter I, is given by the solution of a singular nonlinear 
Volterra integral equation of the second kind. Since the 
integral equation is not in closed form and is nonlinear, 
numerical techniques seem to be the most practical way to 
tackle the problen. Over the past twenty years, a lot of 
research has been done on the numerical solution of an 


integral equation of the form 


U(1,t) = o(t) - Af la + Dy, be FM F[U(1,t)] ae 1) 


Among the existing numerical methods for solving (2.1), the 
method of successive approximations is the most popular one 
(see [1}]). It is based on the idea that the set of successive 


functions defined by 


Your (t) = O(t) - hf "k(t-t) F(¥q(t)) dt, (2.2) 


where k(t-t) is equal to the term in braces in (2.1), 
converges to a solution of (2.1) in every finite interval of 


time. In the following section, the solution method will be 
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outlined, and at the end of the chapter, general comments will 


be made on the technique. 


B. OUTLINE OF THE METHOD 

Consider the time domain in which integral equation (2.1) 
is to be solved. Suppose the domain is partitioned into N 
intervals. For the first time interval, Ost<t, , the 
approximate solution of the integral equation can be obtained 


by using the iteration procedure 


U,.4(1,t) = (t) - hf "k(¢-t) F(U, (1,1) ) dt (2.3) 


until the error between two approximations is less than a 
predefined small number. Next, consider the second time 


interval, t,<tst,. In this interval, 


U..4(1,t) = (t) - hf °k(t-t) F(U,(1,%) ) dt (2.4) 


can be broken into 


U_..(1,t) = o(t) ~ bf °k(t-t) F(U,(1,)) at 
c 
- bf k(t-t) F(U,(1,7)) dt . (2.5) 


Since U(1,t) is determined for O<t<t,, the first integral, 


hf °k(t-t) F(U,(1,4)) at, 
0 
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is known(approximately), and thus the iteration procedure is 
only needed for the second integral. In general, for the i-th 


time interval, the iteration procedure is given by 


U,.(1,t) = $(t) - bf k(e-t) F(U,(1,)) dt 


= hf ° Maee-c)r to (1,t)) av , (2.6) 


where t,.,<tst,. 

As the procedure continues, the surface temperature is 
found for all desired times. One may notice that as the 
algorithm is carried out, the singularity of the Volterra part 
of the integral equation creates difficulty. Appropriately, 
one has to know the nature of the singularity which the kernel 


possesses. To illustrate the idea, consider the integral 


aa dz Cap 


This integral is often found in the integral representation of 
the stated problem. Integral (2.7) possesses a singularity 


which can be removed by the use of the transformation 


Z=at+ (b-a) (1-x7) . (2.8) 
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Then, by using a suitable Gaussian quadrature formula, the 
integral can be evaluated accurately. Normally, one usually 
comes up with an integral with a stronger singularity. 

The iteration procedure outlined above needs a starting 
value. Generally, the algorithm will converge faster to the 
exact solution if the starting value is close to the exact 
solution. Thus, the choice of initial approximation is 
crucial for convergence. Based on the fact that the solution 
of the stated problem is continuous, one can choose the 
temperature at the previous time level as the first 
approximation of the method when a small time step is used. 

The method of successive approximations has been applied 
(in Chapter II) to solve integral equation (2.1). In 
particular, a method used to tackle a simple type of 
Singularity, which one may encounter when evaluating the 
Volterra part numerically, has been discussed. Since numerical 
integration is one of the key steps in the method, the choice 
of the numerical integration scheme does affect the overall 
performance of the algorithm. One can improve the accuracy of 
the successive approximations method by appropriately choosing 
a numerical quadrature that can best deal with the singularity 
found in the integral equation. Even though the procedure 
outlined above may seem simple, it has been shown that the 


method is impractical for large times [3]. 
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IiIt. THE RUNGE-KUTTA METHOD 


A. INTRODUCTION 

This chapter considers another way to deal with the 
integral equation for the stated heat conduction problem, 
namely the Runge-Kutta method which was first introduced by 
Crosbie and Viskanta [5] in 1968. The basic idea of the 
method is based on an approximation of the kernel by a 
separable kernel. The integral equation is differentiated 
with respect to time and transformed into a nonlinear 
differential equation. The Runge-Kutta method is a well known 
numerical scheme for solutions of ordinary differential 
equations. In order to employ the method the surface 
temperature at a desired time must be determined. The order 
of approximation of the method is determined by the order of 
the ordinary differential equation. What differentiates this 
method from the other numerical schemes is instead of solving 
an integral equation directly, the Volterra integral equation 
is first reduced to a system of nonlinear ordinary 
differential equations and then solved numerically. The 
method is not exact since the approximation of the kernel is 
not practical if time steps are small. The accuracy of the 
approximation of the kernel increases with time and order. In 


the next section, the method will be outlined in detail as it 
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is applied to the integral equation (1.72). In addition, as 
an example, the formulas for the third and the fifth order 


versions of the method will be presented explicitly. 


B. OUTLINE OF THE METHOD 

Consider the integral equations derived in Chapter I. 
Generally, the integral representation for the dimensionless 
surface temperature, U(1,t), of the body that we have 


considered can be written as 


U(1,t) = o(t) - f'k(e-t) F(U(2, 7) ) at , (3400) 


where 


k(t-t) = By + WT, Bye ee (3.2) 


The function F(U(1,t)) is the surface heat flux; the a,'s and 
B.'s are eigenvalues and coefficients, respectively. As shown 


in chapter 1, the infinite series k(t-t) has the property 


lim,..K(t) = B, . (3. ay) 
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This is a necessary condition for an integral equation to 
which the method is applied. Now, assume $(t) is a bounded 
differentiable function. The N“-order approximation of k(t-t) 
is given by taking the first N terms of the infinite sum. So, 


(3.2) becomes 


k(t-t) = Bo + Dy Beer” (3.4) 


Substitute (3.4) in (3.1) and let 


I,(t) = et! e**F(U(1, 6) ) dt. (3.5) 


U(1,t) becomes 


U(1,t) = (t) ~ Bof F(U(1,t)) at - Dy, Bele (t) (3.6) 


Differentiating with respect to time, equation (3.6) becomes 


u)(1,¢) =o (t) - B .F(u(1,¢)) 


- ¥, Bgl F(U(1, €)) - wet,(t)) . (3.7) 


yy (2) (1, t) = @ (2) (Cc) _ B Fo) (U(1, t) ) 


ey ees ey) opr, £)) + apz,(t))] , (3.72) 
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u'™) (1, t) = b'™ (¢) = Soe en, cy) 


- DY, BIT, (-1) te Pe (ot, €)) + (-2 "oer, (c) 


In general, the N*-order approximation of the surface 
temperature U(1,t) is determined by assuming that k(t-t) in 
(3.1) takes the form of (3.4), in which only the first N terms 
of the infinite sum are considered, so that U(1,t) is 
approximated by (3.6). Then, by performing N+1 
adifferentiations of (3.6), the resulting system . of 
integrodifferential equations obtained by substitution of 


1, .--, N+tl for m in (3.8) is found to be 
Oe) = Oe) Pro ae) 


=, Pel PAU) eee (309) 


u'2)(1,t) =o) (t) - BF (U(1,t)) 


- >}, 8. LF™ (01, ¢)) - efF(uli,t)) + of7,(t)]) 
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GO (ie) = 6 (t) - oF (u(1, ¢)) 
= eee (-1)fap Fr '2-) (0(1, t)) 


+ (-1)%ef"7,(t)) , (s},alal) 


yy (N+2) (1,t) =  (N+2) (t) - BF” ipihaw t)) 
ete =) oe (U1, t)) 


+ (-1) ao" 7, (oe. (riz) 


To eliminate integrals I,(t), where k=1,...,N, from (3.12), 
consider the first N derivatives of the surface temperature 


which are given by (3.9)-(3.11). Rearrangement of the terms 


in (3.9)-(3.11) yields 
Soo, @ePpty(t) = Oa, £) - 6) (ce) + B,F(U(1,t)) 


+r BF(U(1,¢)) , (3.13) 


- Yr, @tB pT, (t) = 0 (1, ¢) - O(c) + Bor (U(2, t)) 


ype a) seer (oti, c))], (3.14) 
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(-1)" Sy ak’ BT, (t) = 7) (1, £) - o>) (t) + BoP 2 


: pe Bee (-1) 44777 FY 2-4) (U7(1, C) )], (3.25) 


ae 


(-2)¥2 SO ek Ble (the= O™ (1, t) - ce) +pgnl-) (ae 


+ Deer Pid (Heer) (o(1,e))]. (BT) 


In matrix representation we obtain, 


Bay B02 ave B yttin Tea( 6) B, 
-B ya; -B 03 ae —B yoy I, (€) B, 


i) 
— 
WwW 
— 
~] 


(-1) ARG (-1) "B03"? P. (-1) "Ban? Ty., (¢) By-4 
(-1) "1p .02" (-2) B02" .. (-1)%p,02"} | ZO) | | By 


where B,,...,B, are defined as 


Be U8N(1, 6) ag (0b) Bae 


o>, Seu, El (3.18) 
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Be eee Or (tot BF (o(1,t)) 


ao ae Gap om, t)) —agr(u(1,t))) , (3.19) 


oe ert tC) — bc) + Bor?) (u(1,¢) ) 


+ Dy, Brg (72) taki) (oa, £))], (3.20) 


epee (1 Clee aarey + Por 8?) (ui, c) ) 


+ Dyer Boyes (2) FaH ROD (O(a, e))] (3-24) 
Now, let 
Bia; B03 eee B yn 
-B, 0; -B,02 a —B yan 
a ; : eee) 


(—1)ipeage (-1)*pyas... (-1) "Py an, 7 
(-1)" B07" (-1)4*28,a3" we (-1)¥*2B yon | 


By Cramer's rule, I,(t),...,Iy(t) can be expressed as a 


quotient of N x N determinants, given by 
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Iyy(t) = 


B, B02 vee B yon 
B, Baas ” —B yay 


re (- 1) "8 .@2 2N-2 (-1) "Ban? 





_ | By (2) "B03. (-2) BR yal 
Tae Det (A) i 
Bia; B, B303 ae B yon 
“Baas B, -B,03 vse —B yen 


(2) "BaP? By, (2 "Ba? (2) "Bal 





— By (-1) 28,034 (-1) 428 ay 


Det (A) 
Bio; B03 — B, Bay 
Bia -B.02 ws OB, —B yen 


(- 1) "B22 (-1)*Ppi05%? .. By, (-1)"Byen 





(-1) "Bias (-1)"8lo2 ee, (=1) 8B yon 


Det (A) 
2 2 
Bia B02 - By-1%n-a B, 
“Baas -B,03 r -f watts B, 


( 1) *B 42 a? (-1)"B.@>° .. (- 1) "B03! 





(=1)* Bee “i (1) ara oe (51) (enone. 
Det (A) 
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(32233 


(3.24) 


(3.25) 


(3.26) 


where 


Bia; B03 - B yon 

-B,e3 -B03 wee -B Oey 
Det (A) = . : (3.27) 
(1) "Brena 1) eas" .. (=1)*Byan * 


Kal)’ we Boa .. (-1) "Bon 


Next, by using only the fundamental properties of 
determinants, (3.23)-(3.26) can be simplified as a quotient of 


(N-1)x(N-1) determinants (see appendix-A for further detail). 


Define 
Det(A’) = 
- (a@3-@3) - (a3-a?) e — (@y-@3) 
(43-04) (3-03) 2 (ay-a4) 
( -1) N(2N4 -@2""*) ( -1) Ni aon4 -@2""4) & ( -1) N ( ofA -@2""*) 
( -1) N+1 (034 7 -q24"?) ( -1) N+1 (02"? -@2"?) on (-1) N+1 (a 2N-2— 2h?) 


(3.28) 
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The formulas are found to be 


I,(t) = 
(B+03B,) 7 (a@3-«3) aes = ( y-003) 
(B,-a3B,) (a$-a3) :. (ay-@5) 

(By.— (7-1) %a03"*B,)  (-1) 8 (age 4-05" *)  .  (-2) % (aH *- 003" *) 

(By- ( -1) a Salles ( -1) N+1 (as -03"?) ") ( -1) N+1 (20-2 gy 2h2) 





Ba; Det (A’) 


(3.29) 


and for k = 2, ..., Ne-l, 


(continued on next page) 
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ttc) = 


= (a2-aj.;) i - (a2_,-a2.,) (B, +0%.,B,) 
(1 -O fh.) a (h-1 Ops) (B,-a%-1B,) 
(-1)* : sae 
(-1) (apy mage) (28 Cr? -okes) (By. y= (2) Maz, 
(-1) M2 (2h ZRF) (-2 ME ge? 887)  (By- (-1) 40:5 717B, ) 
- (Ohs2-Ahs3) a — (@y-O Ky) 
(Oh. Ops) ae (an-O ke) 


(-1)% (apt Pht) (2) MN 4-254) 
(-1) 992 (p55 837) (2) 88 (7-037) 


Bay Det (A’) 


(3.30) 
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and 


Iy(t) = 


- (ai -a5-1) : — (@y-2-GN-1) (B,+0y-1B,) 


(a3 -Gy-1) A ( y-2-Gy-1) (By -Gy-1B, ) 


(-1) 8 (att -agtt) wn (H 1) (aig -aies*) — (Byy- (-1) Noes" B,) 





(-2) 82 (gy oi?) (2) 2 (ape ates) (By (-1) 82005178, ) 


Bey Det (A’) 
(3.31) 
Thus, I,(t), where k = 1, ..., N, in equation (3.12) can be 
determined explicitly by formulas (3.29) through (3.31). 
Hence, the integrodifferential equation (3.12) is reduced to 


a (N+1)“_ order nonlinear ordinary differential equation 


ot) (ae) = 4) (6) — iB eR ern ore) 


- Dy, Be IS, (-1) tage) (o(1, €)) 


+ (=) *8a3 Sr (t)) , (3.32) 


with the initial values 


U(1,0) = (0) , (3.33) 


u')(1,0) =o (0) - B,F(U(1,0)) 


- }%_, By [F(a 0))) , (3.34) 
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u'23(1,0) = >) (0) - BF (U(1,0) ) 


ie, B,(F (u(1,0)) - agF(U(1,0))) , (3.35) 


u) (1,0) = @™40) — Por (7(1,0) ) 


Py ites re si(1,0))) , (3.36) 


which can be obtained by putting t = 0O in (3.9) through 
(3.12), and I,(t), ..., Iy(t) are determined by formulas (3.29) 
through (3.31). 

The Runge-Kutta method is then applied to the nonlinear 
ordinary differential equation (3.32) with initial conditions 
(3.33) through (3.36). Hence, the N® order approximation of 
the surface temperature will be given by the numerical 
solution of a (N+1)* order nonlinear ordinary differential 
equation. 

As an example, in the next section, formulas for the third 
and the fifth order approximation of the method will be 
presented. Details of the derivation of the equations will 


not be produced, and only major results will be given. 


45 


C. THE THIRD ORDER APPROXIMATION 
The third order approximation of the surface temperature 
of (3.1) is given by the numerical solution of the following 


fourth order nonlinear ordinary differential equation 
yo) (1,t) =o) (¢) - (6, + fehl pe (ae 
+ (Bia + Boa2 + B03) FP) (U(1, t)) 
- (B,a7 + Boaz + B03) PF (U(1, ¢)) 
+ (Piao + B05 + B03) F(U(1, ¢)) 


- (P,ajI,(t) + B,057,(t) + B,057, (Ct) ) (3.37) 


with initial conditions 


U(1,0) = (0) , (3.38a) 


ov) (1,0) =o (0) - (PB, + B, + B. + B,) F(O(0)) , (3.38) 


o'2) (1,0) = @) (0) - (By + B, +B + B,) F' (6(0)) 


+ (Pai + B03 + B03) F(O(0)) , ('383)) 


and 
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oa, O) = gh ao) = (Biome B, + BP, +085) P40) ) 
oh (BP, a3 2 B03 - B03) Fo) (p> (0) ) 


- (B,at + Bja> + Bja$) F(O(0)) . (3.40) 


I,(t), I,(t), and I;(t) in (3.37) are 


peng - (03-03) 
B24 ae 4% 
jE) Se (By762B,) fase) ' (3.41) 


Bia; Det (A’) 





(ai-a3) (B,+a3B,) 


(ai-a$) (B,-a$B,) 


I,(t) = F (3.42) 
; B.a5 Det (A’) 
(a7-03) (B,+03B,) 
4_4 a 
I,(t) =- ea (3.43) 





B03 Det (A’) 


where 


Peer @ cet, + Pp, + PB, + B,)F(u(1,c)) , 
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B, = U2) (17 t) ~My BS B, + be tabs) ae ieee 


+ (Piat + B,03 + B03) F(U(1,t)) , 


B, = 0°)(1,¢) - o°)(e€) + (8B, + 8B, + Bo BS (Cl ieee 
- (B,@; + Boa; + Pye) FUT, £)) 
+ (B,at + Boa> + B05) F(U(1,¢t)) , 


and 


2 2 2 2 
-(@5-a@1) -(a@3-@}) 


DeelA’) = 








alee 4 .4 
(2-3) (03-3) 


D. THE FIFTH ORDER APPROXIMATION 
The fifth order approximation of the surface temperature 
of (3.1) is given by the numerical solution of the following 


sixth order nonlinear differential equation 


48 


atc) — oO (re (6, > 6, +i + B, + B, + 6.) PF ©) (0(1,¢)) 

+ (Bai + Ba2 + B03 + Bag + Bas) FF (U(1, ¢)) 

- (Bier + Boa2 + Bia + Byag + Boas) FS (U(1,¢)) 

+ (B07 + Boaz + Ba; + Bag + Boas) F) (U(1, ¢)) 

- (Bia; + Baz + Ba; + Byag + Bsas) F (U1, ¢)) 

+ (B,a7° + B,a3° + B,a3° + B,ag® + Boas’) F(U(1, ¢) ) 

- (B,@3°2,(t) + B,@3°7,(t) + B,e37Z,(t) 
Pa@e et) + Bas J.(t)) , (3.44) 


with initial conditions 


U(1,0) = (0) , 


u (1,0) = (0) - (By + B, + By + Bs + By + Bs) F(G(0)) , 


u'2) (1,0) = >) (0) - (PB, + B, + B, + Bs + By + BQ) F™ ((0) ) 


A (B04 # B03 + B05 + B04 os: B.as) F((0) ) ' 
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U7) (1,0) = 9) (¢) — (By +B, B+ Bs + Be Be) Ge 
i (Bia; ey B03 a B,03 7 Bia4 + Boas) FF (p (0) ) 
- (Pyar + B03 + B05 + Bag + Bas) F((0)) ' 


yu (1,0) = (0) - (Bo + B, + Bz + Bs + By + Bs) FP ((0)) 


e (B07 ¥ B03 % B05 a B07 a Boag) Fi?) (> (0) ) 


(Bar + B.a2 + B03 + Bias + Boas) F (p(0)) 
+ (Bia: + Boa + Bya3 + Bag + Boas) F(>(0)) , 

u's) (1,0) = (0) - (By + B, + B2 + B3 + By + Bs) FP (h(0)) 

+ (Bar + Boa + B,a3 + Byag + Boas) > (>(0)) 

- (B,a1 + Boa, + Bia3 + Buoy + Boas) F'?) (b(0)) 

+ (Bar + B02 + B03 + Byag + Boas) F ((0)) 


7 (Bia; x B03 v B05 + Bia, Bs Bas) F(p(0)) - 


I,(t), I,(t), I3;(t), I,(t), and I,(t) in 3.44 are 
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2 


(B 


) ( 
5 2 




















ig 
ae) — = (B,-a2B;) (ot3 eta 
37>) (a,- : 
4765) : 
Ba? D aan 
1 Det(A’) : 
| F (3.45) 
- (a7-«2) 
(B,+a5B 
(at-a@3) (B,- hig 
-(a° 3-43B,) (44-03 a 
ao-as3) (By tas al 
reac) (a 4ta53B,) -(ag-a3 sat 
2 = - o-a5) (Be - : ae iw | 
5 “3B,) (as-as ry 
B03 D ey ieee 
5 Det (A’) 
: (3.46) 
-(a3-a3) - (03-02) 
ae ay ‘ (B,+04B,) -(az- 2 
~ (05-28 3-04) «(B,-@4 a 
Gi-G&s) —( 2-04 Fo ne 
eC (a? @2-@,4) (B,+ ; Ae 
ae) = ai -04) 8 ae 
(a5-a2) ( oe 
oat B,-04B, ) (a5—O4 
3%3 Det (A’) — 
4 (3.47) 
-(a7-a2) -(a2-a? 
(a{-as es 
i-@s)  (@3-a5 a 
~(«° o-Ms)  (@3-a5 ae 
af-as) -(aS-ae ct 
Slip || (a? 2-5) -(a°-a' oy 
,(t) = af-as) (as ome: 
man 4 +asB, ) 
5 (B,-a5B,)| 
' (3.48) 


Ba; Det(A’) 


51 


-(@i-a%) - (03-03) -(a3-a2) (B,+07B,) 
(ai-a4)  (@3-a3) (a@3-a5) (B,-@4B,) 


—(a@{-a3) -(a$-a$) -(a3-a3) (B,+a¢B,) 


(as-aa) (O9-G4) (53-04) (B,-agB,) 





5 , (3.49) 
B.a5 Det(A’) 


where 
—(a5-a2) -(a%-a%) -(ag-at) -(as-a7) 
~ (2-01) (3-03) (a g-a) — (@5-a3) 
-(aS-aS) -(aS-0§) -(aS-a8) -(aS-«$)| 
(ao-a?) (aS-a8) (a8-af)  (as-as) 
B, = o™) (1,t) - o) (t) + (PB, + B, + B, + Bs + B, + B.) PCO Re 


B, = u'h(1,t) - o@)(t) + (B, + B, + 82 + B; + B, + Bs) Fie) 
- (Byat + B.03 + B,05 + Bias + Boas) F(U(1,¢)) , 
B, = U)(1,t) - 9?) (t) + (Bo + B, + B. + Bs + B, + B,) F) (U(1, t) 


a (B,a5 = B03 Bs B05 7 B04 + Boas) F‘2) (U( ee 


as (Bia; > B02 v B05 ns Bas % Bas) F(U(1, t)) , 
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eo, Gr etc + (Pp, +p, +P, + 6, + Pp, + PF (O(1, Cc) 
- (Bia; + Braz + Bya3 + Bas + Boas) F? (U(1, t)) 
+ (Byat + B03 + Bja3 + Bas + Boas) FP (U(1,¢)) 
- (Bias + Boa + B03 + Bag + Boas) F(U(1,¢)) , 

and 

B, = US) (1,t) - o>) (t) + (BP, + B, + B. + Bs + B, + BL) F™ (U1, ©) 


- (P,ai + B03 + P5053 + B,0s + Boas) FP (U(1, ¢)) 


~ 


(B,o3 i B03 a B05 i B 0% ‘i Bas) F ') Coa c)”) 
~ (Bias + Boas + B03 + Bag + Boas) FP (U(1, c)) 


(Bia: + Bja2 + Bia3 + Byog + Boas) F(U(1, ¢)) 


+ 


E. REMARKS 

In this work the Runge-Kutta method is applied to solve 
the integral equation resulting from the heat conduction 
problem with combined convection and radiation. In 
particular, the nonlinear ordinary differential equation has 
been determined for both the third and the fifth order 
approximations. It may be observed that this method is not 
very practical for calculating the temperature at small time 


steps. The reason is that the smaller the time one takes, the 
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more terms will be represented correctly, which in turn may 
result in a high-order nonlinear ordinary differential 
equation with a very large number of terms. The number of 
terms could grow to infinity. Thus, the method is usually used 
to compute the surface temperature at large times where the 
temperature distribution is in a steady state. 

In the following chapter, we will take another approach 
using a aifferent numerical method, namely the finite 
aifference method. This method is different from the previous 
numerical techniques in that instead of solving the integral 
equation, it approximates the partial differential equation 


and the boundary conditions directly. 
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Iv. THE FINITE DIFFERENCE METHOD 


A. INTRODUCTION 

The basic idea of the finite difference method is to 
transform a continuous model into a discrete system by 
replacing the continuous domain in the model with a 
denumerable domain. In applying this idea to differential 
equations, all the derivatives in the equation are simply 
replaced by finite difference approximations. Thus, the 
unknowns in the difference equation have a countable domain, 
and the resulting discrete system is solved numerically. 

In the theory of numerical analysis, the significance of 
the computed solution of a finite difference scheme in 
relation to approximating the exact solution depends upon 
three eiements. They are consistency, convergence, and 
stability. Consistency is a condition used to assure that as 
4x (the spacing) approaches zero, the truncation error of the 
scheme also goes to zero. It implies that the finite 
difference can be an arbitrarily accurate approximation to the 
derivative. Convergence of the approximation assures that if 
AX goes to zero, the difference between the computed and the 
exact values also goes to zero. In other words, any desired 
accuracy of the approximated solution can be achieved. The 


last element is the stability. The stability of a scheme 
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concerns the growth of the errors found in the calculations 
which are needed to solve the system of linear equations. A 
scheme is said to be conditionally stable if the roundoff 
error does not amplify if the time step is under a critical 
value which is determined by the differential equation 
considered. In the Lax Equivalence Theorem, the relationships 
of these three conditions are stated. It says that given a 
properly posed initial value problem and a finite difference 
scheme which satisfies the consistency conditions, stability 
is the necessary and sufficient condition for convergence. 
There are many difference approximations and methods for 
solving discrete systems that are available in numerical 
analysis. Different choices of approximation and methods of 
solving the system will lead to differing degrees of accuracy 
in the approximation of the solution. This chapter will only 
focus on a particular finite difference scheme used to 
approximate the governing partial differential equation in the 
stated problem and an algorithm for solving the discretized 


systen. 


B. CRANK-NICHOLSON SCHEME 


Suppose a lies between x, and x; and t 2 t,, where x, and 


t 


xX; are some initial and final x-coordinate which brackets the — 


location of concern. Let ax and at be increments of x and t, 
respectively. The x-t space can be partitioned into a grid 


network in which the points are given by x = x, + jax and t = 
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t, + nat, where j = 0,1,2,...,N, with N being the number of 
nodes, and n= 0,1,2 ... . When ax and at are constants, the 
mesh obtained is uniform, and the temperature at x = x, + 
iax, written as x,, and t = t) + nax, written as t,, is denoted 
by U,". 

As previously mentioned, there are several ways of 
choosing a finite difference operator for replacing the 
derivatives. If the average of the forward and backward 
difference schemes is used for the space discretization and 
the forward difference scheme is written about the point x;, 
t,.,, the governing partial differential equation becomes a 
second order accurate (in both x and t) finite difference 


equation. It is given by 





Up, + (-2 - 2B) U7" + uf, = - Uf, + (2 - 2B)0P - of, , (4.1) 
where 
= Ax? 
B at ‘ 


(which is the well known Crank=-Nicholson scheme). 

Since it is of second order, the truncation error 
associated with (4.1) is on the order of o(4x’ + at’). Notice 
that the temperature at time t,,, is a function of unknown and 
known temperatures at six of the ten points shown on the Fig. 


4.1. 


Si 





i-l 1 itl 


Fig. 4.1 Graphical representation 
for Crank-Nicholson Equation 


Though the Crank-Nicholson scheme is’ unconditionally 
stable in the von Neumann sense [10], it is well known that 
the finite difference solution may suffer from severe 
oscillations when it is applied to the Robin's type boundary 
conditions with the choice of a time step which is large 
relative to the spatial step [10]. This phenomenon becomes 
more pronounced especially when the nonlinear boundary 
conditions are treated with a crude approximation where the 
surface temperatures are not determined with sufficient 
accuracy. Thus, in the problems considered, the Crank- 
Nicholson scheme is only applied to the boundary at x = 0 and 


the space between the boundaries, that is , 
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OG, =e, N-ae for case 1 
if | = 


1, ..., N-1 for case 2 


To ensure that the oscillation is eliminated, the implicit 
backward finite difference scheme (which is satisfactory with 
all types of boundary conditions) is adopted at the boundary, 


Xx = 1. The equation at x = 1 is given by 


Ue. + (2 - PU + Ue =~ BU, , (4.3) 


where $ is as before. 

There is a fictitious point outside the computational domain 
in (4.3), that is, the unknown temperature at N+1 is denoted 
as U,,,""'. To eliminate that point, use a difference method to 
approximate the derivative in the radiative boundary condition 


(1.2c) because 


n+1 n+1 


Une1 ~ Uy- + + 
Ti TER. eyog =r(ug) (4.4) 


where F is the right hand side of (1.2c). Algebraically 


manipulating (4.4) yields the following equation 


User = Uner + 2ax0,Uy + 2axF (Uy) . (4.5) 


Substituting (4.5) into (4.3), the resulting expression 


becomes 
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2Uy-y + (-2 - B + 2axa,) Uy” = - BUYS - 2axF(UN"*) , (4.6) 


which is a nonlinear equation in U,"*:. 

Observe that (4.1) and (4.6) constitute a set of 
simultaneous equations at each time. step. In matrix 
representation, the resulting system is of the form 

AU = B : (4.7) 
where A is a tridiagonal matrix, B is a vector of all the 
known values found in each equation, and U is a vector of the 
unknown temperatures at each space node at a particular moment 
of time. So, for each time level, the transient temperature 
is given by the solution of a system of equations. 

The Thomas algorithm can be used to solve a tridiagonal 
system of linear equations. Clearly, all the equations in 
(4.7) are linear (except the last one). The first half of the 
algorithm, as given in appendix-B, can be directly applied to 
the system except for the case where i=N. In the case of i 
= N, substituting d(N) in the first Do-loop in the expression 


right after the first loop yields 


d(N) - ratio * d(N-1) (4.8) 


n+1i _ 
Un b(N) 
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which implies 


-BU,y - 2axF(Uy**) - ratio * d(N-1) 


yar = 
sa b(N) 


(4.9) 


with b(N), d(N-1), and ratio computed in the first Do- 


loop(reference to appendix-B). Now, (4.9) can be rewritten as 


-BU, - 2axF(Us"*) - ratio * d(N-1) 


_ ppntl _ 
Sie we 0. (4.10) 


Let the left hand side of (4.10) be represented by f. It 


follows that 


AiG )r= 0 (An) 


Thus, the update of the surface temperature is the solution of 
the nonlinear equation (4.11). 

In the following section, the cases for a flat plate and 
a sphere will be considered to obtain the respective 


tridiagonal systems. 


C. TWO SPECIAL CASES 
1. The Flat Plate 
The parameters corresponding to this case can be found 


in Chapter I. Applying the finite difference method outlined 
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above to the governing equations leads to the following 
results. 


Consider the Crank-Nicholson scheme for j = 0, ..., N-1. 
For j = O, 


UN + (-2 - 2B) ug + uf = - U0 + (2-2B)Uu27 - U2. (4.12) 


To eliminate the fictitious points, the boundary condition at 


xX = O in discretized form is taken to be 


n+1 n+1 
U; aa UL; 


= : 4. 
oe 0 (4.13) 
Thus, 
uy’ = Un” (4.14) 
and 
UY = U4 (25 ee 
‘ 
| 
Substituting (4.14) and (4.15) in (4.12) produces 
(-2 - 2B) US" + 2u7"* = (2-28) Ug - 2uU2 . (4.16) 
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Bote =p ..., N=1l 


Uy + (-2 - 2B) Up" + us" = - uf + (2-28) U2 - US, 


Uy-2 + (-2 - 2B) Uy + Uy = - Uyn + (2-28) Uy, - UP . 


When j = N, as shown before, the equation becomes 


2Uy-, + (-2 - B + 2axa,) UN" = - BUy - 2axF(UX"*) 


In matrix representation, with initial values 


Uj; =1, where j=0,...,N, 
we have 
(-2-26) 2 0 0 
zl (-2-26) 1 0 0 
0 il (-2-28) 1 0 
2 0 0 < 0 
a e e 0 
0 0 eee 1 {(-2626) 1 
0 0 eee 0 2 (-2-B+2axa,) 
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a 


ee i) 


(4.18) 


(4.19) 


(2-25) Uo. = 208 ae 


-U> +2528) OF - we Utne 

hs | , and u=| |. (4.20) 
-Uy-2 + (2-2B) Uy. - Uy Un-1 
-BUy - 2axF(Uy"*) On 


2. The Sphere 
The parameters again can be obtained in Chapter I and 
will not be repeated here. Using the finite difference method 
outlined above with the governing equations leads to the 
following: 


Consider the Crank-Nicholson scheme for j = 1, ..., N-1l. 


For j = 1, 
Us * + (-2 - 2B) up** + us" = - uy + (2-2B)Uf - UP. (4.21) 
However, 
Us * =0 and uU;=0 . (4.22) 


Substituting (4.22) in (4.21) produces 


(-2 - 28) uf + uy" = (2-2B)UZ - Uy. (4.23) 
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C2 — 20)0, + U, = -9U, + (2=28) uy - vu; , (4.24) 


Ung + (-2 - 2B) Uyey + Uy™* = - Uye + (2-28) Uy - Uy . (4.25) 
When j = N, as shown before, the backward scheme is used. 
Thus, 

2Ug) + (-2 - B + 2axa,) UV" = - BUY - 2axF(Uy) 


In matrix representation, with initial values 


U; = jax, where j=1,...,N, 
we have 
(-2-26) at 0 0 
al ( <2e-2(3) ad: 0 0 
0 1 (2-2 0 
0 0 0 
A= . 5. ata, 26) 
: : 0 
0 0 i {=2=26) il 
0 0 0 2 (-2-B +2axa,) 
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(2-28) Te _ ie of 


=U; + (@G20ii 0, ne 

Jz = } , and u=| | |. im 
-Uy-2 + (2-28) Uy. - Uw Un-1 
-BUS - 2axF(Uy"*) Ux 


D. STABILITY 

Even though the backward analog is implemented on the 
radiative boundary, according to the numerical experiments, 
the method still suffers from the problem of oscillations when 
a large time step is imposed. As far as the author is aware, 
not a single formula has been developed for the stability 
criteria of an implicit finite difference scheme with 
nonlinear boundary conditions. However, two stability 
formulas of the related problems, which are found in the 
literature [7], can serve as a guideline in choosing the time 
step for the problems considered. The first one is due to 
Lawson and Morris [7]. They deduce the stability criterion 


for the Crank-Nicholson equation with linear boundary 


conditions as 


at < a (4.28) 


Another stability criterion is due to Milton and Goss [9] who 


applied the laws of thermodynamics in developing the stability 
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requirement for an explicit finite difference scheme with 


nonlinear boundary conditions. It turned out that the time 


step required for the stability is restricted by 


he (ax) ?max{u,7} 


a (4.29) 
Beee> Ue = bX(Uj)* - «,axUy} 


where the maximum is taken over all n and where U,° can be 


found by setting the following function 


AU 
ae = A — BU, 


at CCU)" 


equal to zero, where 


— pypotl _ n 
40, a On On 


f 


_ 2Uyey 


(ax) 2 


— lax = 1) 
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Because this formula is not very practical in actual use, 
(4.28) will be chosen as a guideline for selecting the time 


step of the method. 


E. REFINEMENT OF PARTITION AND EXTRAPOLATION TECHNIQUES 

The partition of the domain covered has a great influence 
on the accuracy of the solution obtained. The choice of grid 
points is determined by knowledge of the problem and by 
numerical experimentation. Here, two ways of improving the 
accuracy of the finite difference method are presented, and 
one of the two is chosen to be implemented in the numerical 
methods. 

One way to improve the accuracy is the so called 
prolongation. We first solve the problem using one spacing 
and then refine the partition and then repeat the computation. 
If the comparison shows large differences, the process is 
repeated for smaller and smaller grid sizes until a desired 
accuracy is achieved. This method may result in a prolonged 
computational time for the solution. 

The second way is called extrapolation. The simple 
ingenious idea of the technique, which dates back to 
Richardson in 1910, is the following: 

One solves the same type of problem over a prescribed 
interval, for example [0,1], several times with successively 
smaller step sizes. Thus, one obtains a sequence of 


approximations 
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ace, Y\lyia,) ; 


for a given sequence of step sizes 


ee Os 


The successive step size h, is often defined in terms of an 


input step size h by 


h;=—, 10,1,2,...  . (4.31) 


Thus, any step size sequence {h,} can be characterized by the 
associated integer sequence {n,}. The following are some 
examples of integer sequences: 
mar, 2a, SHG Sa, ...)} (Romberg sequence) 
(2 246 , 812 2s} (Bulirsch sequence) 
{1,2—604,...)} (harmonic sequence) 
So, the numerical solution at x is computed for a sequence of 
step size h, and denoted by 1T,, = y(1,h,;).- Then, the 
extrapolation tableau, 
Too 
Tio Tis 
T20/T21/T22 
is calculated for x according to two types of commonly used 
extrapolation schemes 


a). Aitken-Neville algorithm 
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For i = 1;2,... anaek E2737 


= Ts ,-1 = Ty, ke1 ai P3-1,k-1 





i,k : (4.32) 
oe 3 
b). Rotational extrapolation 
For i =oy2 ,...@and k. = 92, 3y4aee , 2 
Ti, = 9 ' 
a = 1 jae 
Aes ie faa + 1,K-1 1-1,k-1 (4 33) 





2 
n; —_ Dy,K-a 7 Ti-1,K-1 = 
Dy, k-1 _ Dy-1, k-2 


In this study, extrapolation scheme (4.32) with the Romberg 
sequence will be implemented when the finite difference method 
is RSeaKes find the numerical solution of the stated problems. 
It should be noted that if this extrapolation scheme is used 


the computational time will increase exorbitantly. 
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V. NUMERICAL RESULTS 


A. INTRODUCTION 

The problem described in section 1(B) was’7 solved 
numerically for two special cases, namely, the flat plate 
(@,=1, @=0, a@,;=-1, h=1) and the sphere (@,=0, @,=-1, @,=1, 
h=1). Since a lot of numerical results of the problem 
computed by successive approximations method are available in 
some of the papers[(3,4], in this thesis, only the Runge-Kutta 
method and the finite difference method are employed to the 
problem for study. Programs are written in Fortran 77 using 
the Amdahl 5990 model 500 mainframe computer and are set up-to 
allow input for the time step. Thus one can approximate the 
maximum time step that can be used in a particular numerical 
method. All calculations are done using double precision 
arithmetic yielding 12-digit accuracy. Numerical results 
generated by the methods are compared and discussed. 

The Runge-Kutta and the finite difference methods are 
implemented to solve both special cases. In particular, three 
different order approximations of the Runge-Kutta method are 
programmed to solve the integral equations derived by the 
Laplace transform method. Inefficiency of a high order Runge- 
Kutta method motivates the use of the finite difference 


technique. Again, the method is implemented in both cases for 
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various time steps. Some of the numerical results are 
tabulated and plotted in such a way that a comparison can be 
made. Notice that the Runge-Kutta method is not applied to 
the integral equations obtained by the eigenvalue expansion 
method. The reason is the lag parts of those integral 
equations diverge when time is zero, and thus, the initial 
values of the nonlinear ordinary differential equations cannot 


be computed. 


B. RESULTS FOR THE FLAT PLATE AND THE SPHERE 

Integral equations (1.39) and (1.43) are solved using the 
Runge-Kutta method of orders 1, 3, and 5. The first order 
approximation can be found in [5], whereas the third and the 
fifth order approximations are described in sections 3(C) and 
3(D), respectively. Solutions of the nonlinear ordinary 
differential equations corresponding to (1.39) and (1.43) are 
obtained using the fourth-order Runge-Kutta method developed 
by Zurmuhl] [15]. The results show that solutions of a high 
order approximation fall below those of a lower order 


approximation (Fig. 5.1, Fig. 5.2, Fig 5.3). 
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0.9 —-— ist order 
0.8 \ — — 3rd order 
0.7 
0.6 
0.5 
0.4 





time 


0.2 0.4 0.6 0.8 i. 


Fig.5.1 Surface temperature of a Fiat Plate 
cooled by convection and radiation 
(Runge-Kutte Method). 
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Fig.5.2Surface temperature of a Flat Plate 
cooled by convection and radiation 













(Runge-Kutta Method). 
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Fig.5.3 Surface temperature of a sphere 
cooled by convection and radiation 
(Runge-Kutta Method). 
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With respect to time step, we do not have the same phenomenon 
as in the order of approximation. In a fixed order 
approximation method, the solution curves for a smaller time 
step fall below those for a larger time step at small times 


(approx. less than 0.2) and above at large times (Fig. 5.4). 


Surface Temp 
0. 91\ —— 0.01 
0.8} \ = ac.) 
0.7 
0.6 
0.5 
0.4 





Fig. 5.4Surface temperature of a Flat Plate 
cooled bY convection and radiation 
(Runge-Kutta Method of the First Order) 


According to numerical experiments, the stability requirements 
for Runge-Kutta method of orders 1, 3, and 5 are approximately 
0.1, 0.01, and 0.001, respectively. As observed earlier, a 
drawback of the Runge-Kutta method is that it requires the 
solution of a heuristic nonlinear ordinary differential 
equation for a high order approximation. This leads to an 
attempt to use an easier algorithm, and for this reason, the 
finite difference method was implemented. 

Equations (1.1) and (1.2a,b,c) are solved for the flat 


plate and for the sphere. The extrapolation formula used to 
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improve the accuracy of the solutions is the Aitken-Neville 
algorithm (4.32). The results for various time steps are 
presented in tables 1 and 2, and some of these results are 
plotted in Figures 5.5 and 5.6. As tables 1 and 2 show, the 
situation where solution curves for a smaller time step fall 


below those for a larger time step holds in the finite 


aifference method. 


Table i 
The Finite Difference Method for Various Time Steps 




















































































0.01 0.849395 0.843059 0.842539 
0.02 0.797214 0.793829 0.793609 
0.03 0.764701 0.762365 | 9.762256 
0.04 0.740419 0.738609 0.738559 
0.05 0.720827 0.719340 0.719328 
0.06 0.704316 0.703048 0.703063 
0.07 0.689998 0.688893 | 0.688927 
0.08 0.677333 0.676351 0.676400 
0.09 0.659595 | 0.665075 | 0.665137 
0.10 0.655626 0.654821 0.654893 
0.20 0.583564 0583127 
0.30 0.535891 0.535564 
0.40 0.496573 0.496292 
0.50 0.461310 0.461056 
0.60 0.428810 0.428575 
0.70 0.398614 0.398394 
0.80 0.37050! 0.370296 

0.344321 0.344130 





0.319948 0.319770 


The Surface Temperature of a Flat Piate cooled by 
Convection and Radiation. 
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0.6 
0.4 Seliatee Sak 
0.2 
; time 


0.2 0.4 0.6 0.8 1. 
Fig.5.5 Surface temperature of a flat plate 
cooled by convection and radiation 


(Finite Difference Method, At=0.0 4) 


Tablé 2 
The Finite Difference Method for Various Time Steps 















0.916249 | 0.912694 0.913001 




































0.02 0.882656 | 0.880526 
~ 0.03 0.860270 | 0.858708 

0.04 0.842779 | 0.641514 

0.05 0.828134 | 0.827055 

0.06 0.815380 | 0.814431 

0.07 0.803988 | 0.803134 

0.08 0.793626 | 0.792846 

0.09 0.784075 | 0.774505 

0.10 0.775178 

0.20 0.706513 

0.30 0.656561 

0.40 0.616906 

050 0.584412 

0.60 0.557215 

0.70 0.534052 

0.80 0.514033 

0.90 0.496514 


0.481017 


Thé Surface Temperature of A tphereé cooled by 
Convection and Radiation. 
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Surface Temp 
0.02 0.04 0.06 0.08 
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0.1 
0.95 
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0.9 a — O04 


0.85 


0.8 


Fig.5.6Surface temperature of a sphere 
cooled by convection and radiation 


(Finite Difference Method). | 


Even though the implicit scheme is implemented on the 
boundary, numerical experiments show that solutions still 
exhibit oscillation when a large time step was chosen (time 
step > 0.01). This constraint of time step leads to large 
computational times for large time solutions. 

Figures 5.7 and 5.8 show representative results for the 
Runge-Kutta and finite difference methods where a flat plate 
and sphere are cooling. Tables 3 and 4 show that, when at = 
0.01, the results obtained by using the Runge-Kutta method of 
orders 2 and 3 compared favourably with those using the finite 
difference method. The difference of the solutions by using 
the two methods is less than 3.1% (relative error) in average 


for each case. 
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Table3 


Comparison of the Runge-Kutta and the 
Finite Difference Methods (At=0.01) 


Time | ist Order | 3rd Order] Finite 


Diff. 
0.10 0.712925 0.667017 0.655626 
020 0.619477 0.567586 0.583564 
0.30 0.562454 0.514915 0.535891 
0.40 0.516209 0.477480 0.496573 
050 0.475091 0.445843 0.461310 
0.60 0.437583 0.417183 0.428810 
0.70 0.403124 0.390595 0.398614 
0.80 0.371405 0.365742 0.370501 
0.90 0.342189 0.342459 0.344321 
1.00 0.315274 0.320636 0.319948 


The Surface Temperature of a Flat Pilate cooled by 
Convection and Radiation. 


Surface Temp 


—-— Istorder 
0.9 —-—— 3rdorder 
0.8 —— Finite Diff. 





time 


0.2 0.4 0.6 0.8 i 
Fig.5.7Compersion of results for cooling 
@ flat plate. 
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Table 4 


Comparison of the Runge-Kutta and the 


Finite Difference Methods 


Time 


0.10 
0.20 
0.30 
0.40 
0.50 
0.60 
0.70 
0.80 
0.90 
1.00 


ist Order 


0.785102 
0.705908 
0.651471 
0.609637 
0.576138 
0.548557 
0.525344 
0.505453 
0.488156 
0.472926 


3rd Order 


0.754306 
0.665482 
0.622430 
0.591087 
0.565188 
0.542988 
0.523640 
0.506582 
0.491399 
0.477772 


Finite 
Diff. 


0.775178 
0.706513 
0.656562 
0.616906 
0.584412 
0.557215 
0.534052 
0.514033 
0.496514 
0.481017 


(at=0.01) 


The Surface Temperature of aSphere cooled by 
Cortvection and Radiation. 


Surface Temp 









time 


———e 





0.2 0.4 0.6 O.8 le. 
0.9 
o st * "wise Order 
\ ———- 3rd order 
0.7 —— Finite Diff. 
0.6 
0.5 
(At=0.01 ) 
0.4 
Fig. 5.8 Comparsion of results for cooling 
of asphere. 
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Lastly, in Figures 5.9 and 5.10, the results of two special 
cases solved by finite difference method are compared. The 
graph shows that, when At = 0.01, the surface temperature of 
a flat plate fell much faster than that of a sphere. We 
believe that it is due to the effect of the boundary condition 
at x = 0O and the difference in the coefficient of the 


convective term. 


Surface Temp 
0.3 
0.25 
0.2 
0.15 
0.1 
0.05 


en” er ee en, 


Fig.5.9Surface temperature of a flat plate 
cooled by convection and radiation 


(Finite Difference Method, At=0.0 1) 
Surface Temp 
0.8 
0.6 


0.4 


0.2 
2 ¢ 6 = 3 0 uae 


Fig. 5.10Surface temperature of a sphere 
cooled by convection end radiation 


(Finite Difference Method. At=0.0 1) 
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VI. CONCLUSIONS 
The study of the one dimensional heat equation subject to 
combined convective and radiative boundary conditions in 
rectangular coordinates is motivated by the advent of space 
technology where knowledge of the temperature of bodies in 
deep space is necessary, for instance, in the design of space 
shuttles. 

The solids are assumed to be homogeneous, isotropic, and 
opaque to thermal radiation and to have temperature 
independent physical properties. This assumption leads to a 
linear heat equation. The difficulty of the problem is 
determined by the conditions prescribed at the Boundaries. 
According to the laws of physics, the heat flux of the 
radiative heat transfer is proportional to the fourth power of 
the temperature which causes nonlinearity at the boundaries. 

Problems of this type are first solved by analytic 
techniques, one of which is the integral transform method. In 
particular, Laplace transform and eigenvalue expansion are 
used. The solutions which are explicitly determined at the 
surface for two special cases, namely, the flat plate and the 
sphere, are singular nonlinear Volterra integral equations of 
the second kind. Although they are not practical in 


determining the temperature at a particular time, these 
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integral equations can help us to deduce some useful 
information about the behaviour of the surface temperature. 

Since the analytic solutions found for the problem are not 
practical to use, numerical techniques are considered as an 
alternative. Two numerical schemes that are used to deal with 
the resulting integral equation are the Runge-Kutta method and 
the successive approximations method. Both techniques are 
studied in great detail. Numerical solutions show that the 
successive approximations method is "exact" in the sense that 
any desired accuracy may be obtained [3,4]. Additionally, the 
closer the initial approximation was to the exact solution, 
the faster the method of successive approximation converged to 
the exact solution. Conditions for the numerical solution and 
limitations of these schemes are also discussed. 

Another numerical technique which is directly applied to 
the governing equations is presented as a possible alternative 
to the numerical methods previously discussed. It is the well 
known finite difference method in which the Crank- Nicholson 
scheme, the backward implicit scheme, and the Newton-Raphson 
method are combined to solve for the surface temperature. 

The Runge-Kutta methods of orders 1, 3, and 5 are 
programmed for (1.39) and (1.43) which are the integral 
equations corresponding to the flat plate and the sphere, 
respectively. The numerical results are presented with 
respect to their orders and to their time steps. The data 


reveal the following phenomena. First, the solutions of a 
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high order approximation fall below those of a lower order 
approximation. This phenomena is a result of the higher order 
approximations closing in on the actual solution. Second, the 
first phenomenon does not occur in the solutions for various 
time steps with a fixed approximation order. The main result 
here is that a smaller step size determines the surface 
temperature for very small times (0<t<0.2) more accurately and 
a larger step size determines the surface temperature for 
larger times (t20.2) more accurately. Third, the agreement 
between the 12%, 322, and 5 order Runge-Kutta approximations 
is better for the sphere than that for the plate. Physically, 
this is due to the fact that the boundary surface area to 
total volume ratio is largest for the sphere and smallest for 
the plate. The reason for this trend is that the larger the 
ratio the more uniform will the temperature be throughout the 
body. fhe fourth phenomenon is that the accuracy of the 
approximation increases with time. For large values of time, 
the rate of change of temperature is reduced, as would be 
expected from the influence of the fourth power term (U’). 
Since the Runge-Kutta method did not offer any efficiency in 
the area of high order approximations, the finite difference 
method was considered. 

Equations (1.1) and (1.2a,b,c) are solved numerically 
using the finite difference method for both the flat plate and 
the sphere. The results for various time steps are presented. 


The table shows that the second phenomenon found in the Runge- 
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Kutta method again occurs in the solutions generated by the 
finite difference method with respect to time step. Again, as 
in the Runge-Kutta method the smaller step size determines the 
surface temperature more accurately for small times and the 
larger step size determines the surface temperature more 
accurately for large times. 

Finally, two comparisons are made of the numerical 
solutions. The first is of the Runge-Kutta method and the 
finite difference method. The results show that there is a 
good agreement between the two methods, and the difference 
between their solutions are, on the average, less than 3.1% in 
both cases. The second comparison was made between the 
solution of a flat plate and that of a sphere. The finite 
difference method conveys that temperature of a flat plate 
decays much faster than that of a sphere. This result was 
expected for the transient heat conduction with linear 
boundary conditions. This could be due to a larger area on 
the plate exposed to the uniform’ boundary layer. 
Additionally, this result could be caused by the sphere having 
a larger surface area to volume ratio; thus the sphere would 
have a more uniform temperature distribution throughout the 
body resulting in a slower decay of surface temperature. 

Comments of a more general nature are included. 

1. The convection mode of heat transfer appears to be 
dominant as the dimensionless temperature approaches 


uniformity for a plate cooling to a zero environment. This 
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result is due to the fact that U* is approaching zero at much 
a faster rate then U. 

2 Physically, the adiabatic or initial temperature 
cannot be equal to absolute zero, however, in many situations 
the temperature ratio of adiabatic surface temperature to 
initial temperature can be very small. 

3. For cooling and heating the solutions are initially 
inaccurate due to the fact that at t=0 the linearized heat 
flux is not equal to the actual flux. 

4. For a set time step size the number of iterations 
required to meet a set accuracy is determined by which 
surface is receiving the highest heat rate. 

5% The time required to achieve a particular surface 
temperature during cooling decreases as the ratio of the 
environment temperature to initial solid temperature 
increases. 

To conclude this thesis, a numerical scheme is proposed as 
an alternative to the existing numerical methods. The method 
of successive approximations is described in Chapter II. One 
of the major difficulties of that method is choosing the 
initial approximation for the iteration procedure. As 
mentioned earlier, the convergence of the algorithm can be 
accelerated if one could obtain an initial approximation which 
is close to the exact solution. To determine this value, one 
could first use the finite difference method (without the 


extrapolation algorithm) described in Chapter IV to determine 
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the surface temperature. Then, by treating it as an initial 
approximation, the method of successive approximations is 
applied to obtain the _ solution. We believe that the 
temperature obtained by using the finite difference 
approximation for the exact solution would be a better 
solution than the temperature at the previous time level. In 
addition, this technique would allow larger time steps. 
However, as far as the author is aware, nothing has been 
proved for this method, and the analytical and numerical 


justifications for the algorithm are left open. 
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APPENDIZ<-A 


To provide a better understanding of the results for I,(t) and 
I,(t) corresponding to (3.29) and (3.30) respectively, 


consider the following example where N = 3 in (3.17): 


Bia; B.03 B05 u. (fC) B, 
-B,03 —B,02 —B305 T,(t)| = |B (Aga) 


Bias B05 B05 I, (t) ae 
L J 


where 


Bia: B02 B03 
A =|-B,a: -B,02 -B,a3| - (A.2) 


Bias B03 B05 


Using Cramer's rule 





I,(t) = (A.3) 


Det (A) 


where 
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Bia; B02 B03 
Det (A) =|-B,a; -B,a2 -B,a5] . (A. 4) 


Bias B03 B05 


Using a fundamental property of determinants (A.4) can be 


written as 


oe 
2 2 2 
Det(A) = -B, BB, (afaza>)|%1 %2 3 . (A.5) 
Gi a2 a3 
Column reduction gives 
1 0 0 
2 2 2 2 2 
Det(A) = -B,B,B,(aiaza5)\%2 ~(%2-a3) —(@3-@3)] . (A.6) 


a 4 4 4 4 
@; (@2-a@,) (a@3-a;) 


Define 
1 0 0 
Al =|az -(@2-a7) -(a3-at)] . wars 
a: (@2-@{)  (a$-a}) 
then 
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-(@3-a3) -(a3-a}) 
Det(A’) = a a a. (Aye) 
(02-G3) (03-3) 
By a Similar manipulation (A.3) becomes 
(B,+03B,) -(a3-03) 
4 4 4 
B,-45B 3-a 
Te C) = { = ) (43 ~ G2) (A.9) 


Bia> Det (A’) 


After working through a considerable amount of algebra both 


(A.3) and (A.9) give 


B,aza2 + B(as + as) +B 
7G) 142-3 2 3 2 3 


27,2,2 2 2 
Biaerleze3 - ai(a3 + a3) + a4] 


In a Similar manner I,(t) and I,;(t) can be determined. Without 


loss of generality I,(t) can also be found. 
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APPENDIZ=B 


The Thomas Algorithm 


The equations are: 
+ C,U;,, = dy , 


a,U;., + DU; 
Oo, and N is the number of nodes 


where 1 < i < N with a, = Cc, = 


in the domain. 


The algorithm is as follows: 


2,N 
a;/b, 

b, =i ratio * Cy-4 
d, a ratio * q,-1 


Ky 

pel] 

ct 

fe 

O K- 
it il not 


10 CONTINUE 
Uy = Ay/by 


DO 20 i N= 
(dy - Cy * Gy) /B; 


Uy 


20 CONTINUE 
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